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Abstract 

We define a probabilistic scheme to compute the distributions of periods, tran- 
sients and weigths of attraction basins in Kauffman networks. These quantities are 
obtained in the framework of the annealed approximation, first introduced by Derrida 
and Pomeau. Numerical results are in good agreement with the computed values of 
the exponents of average periods, but show also some interesting features which can 
not be explained whithin the annealed approximation. 
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1 Introduction 



Kauffman Networks are a disordered dynamical system. They were introduced in 1969 as 
a simplified model of the genetic regulatory system acting in cell differentiation |]]. The 
debate about cellular automata grew the interest of physicists about the model @. It 
was soon recognized that it has deep analogies with a disordered system studied in 
statistical mechanics, the infinite range Spin Glass, whose properties had been investigated 
and understood in the immediately former years ||. 

A configuration of the system consists of N binary variables, (7j = 0, 1, % = 1 • • • N. In 
the biological metaphor, the variable <7j represents the state of activation of the i-th gene in 
the cell. 

The configuration space is provided of a metric structure through the Hamming distance 
(normalized to one), say the fraction of variables whose states are different between two 
configurations C and C: 

1 N / n — n' \ 2 

d(C,C') = ±J2(^i) . (1) 

The system is disordered in the sense that the dynamic is deterministic, but its rules 
are chosen at random at the beginning in the following way: for each gene i, K genes are 
extracted randomly to send a signal to it (the same gene can be extracted more than once); 
these are denoted by I = 1 • • • K; for each gene i and each of the 2 K configurations of 
the signal it can receive, the value of the response function fi{a\ ■ ■ ■ ax) is extracted to be 
with probability p and 1 with probability 1 — p. The dynamic rules are then quenched and 
do not change during the evolution of the system. 

At each time step the N variables are updated simultaneously, according to the signal 
received and to the response to that signal: 

Gi{t + l)=fi (<7 il(i ), • ■ ■ (Tj K (i)) . (2) 

In a more compact way, the evolution law can be put in the form C(t + 1) = $> v (C(t)), 
where C denotes a configuration of the dynamic variables and rj a configuration of the 
dynamic rules. 

Due to the fact that configuration space is finite, every trajectory ends up, after a tran- 
sient time, on a periodic orbit, and configuration space is partitioned into a certain number of 
attraction basins of weight W a , where the a labels the periodic orbit. So, as it is customary 
in the theory of disordered systems, three kinds of average have to be defined: 

• The time average, which is equivalent to the average over a given limit cycle. 

• The average over the same network, where different cycles have to be weighted with 
their own W a . 

• The average over different realizations of the network. 
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It was soon discovered that Kauffman Networks have two different kinds of dynamical 
behaviours: a frozen regime, when the average cycle length grows as a power of N, and a 
chaotic regime when it grows exponentially with N. 

It is possible to characterize this transition as a phase transition and two order parameters 
have been found to describe it. 

Derrida and Pomeau considered the average limit distance between configurations on 
two randomly chosen trajectories in the same network. They computed its value in the 
framework of the annealed approximation and found that it is zero in the frozen phase (in 
the limit of an infinite network: remember that the distance is defined as the fraction of 
different variables), while it is different from zero in the chaotic one [§], || [7|. 

Then Flyvbjerg introduced the concept of stable core, defined as the set of the variables 
that evolve to a constant state not depending on the initial configuration. The fraction of 
variables belonging to the stable core is 1 (in the limit of infinite N) in the frozen phase 
while is less than 1 in the chaotic phase [[13] . 

The transitions described by these parameters take place at the same critical point, 
p = p c {K), where p c is such that 2p c (l - p c ) = l/K g, §, 0, fTJ. 

Despite of the success of this approach, the connection between these order parameters 
and dynamical properties (of the kind of the average cycle length or the weights of attraction 
basins) is still lacking. This work is a contribute to fill in this gap. 

In the second section we define the closing probabilities, which are the quantities that 
allow to compute the full cycle length and transient time distribution. Section 3 describes 
the annealed approximation, first introduced by Derrida and Pomeau M, that here is used 
to compute the distribution of the overlap between different time configurations and the 
closing probabilities, whence we obtain the distribuition of cycle lengths (section 4) and 
the distribution of attraction basins weigths in the chaotic phase (section 5), which turn 
out to be the same as those obtained by Derrida and Flyvbjerg in the limit case K — > oo 
(the so called Random Map) ||. Section 6 deals with the average number of cycles in a 
given network, which is computed in the chaotic phase, in the framework of the annealed 
approximation. 

In section 7 numerical results are presented. This section is subdivided into three sub- 
sections: in the first one, data concerning the distribution of the overlap are shown; they 
support the validity of the annealed approximation, but only in some range of the temporal 
variables. In section 7.2 simulations about the closing probabilities are reported, which show 
deviations from the annealed approximation after an initial agreement. In section 7.3 we 
show the exponents which give the N behaviour of the average cycle length for different 
values of the parameter K; there is a good agreement between our computations and the 
numerical results, but the distribution of cycle lengths that we measured is quite different 
from the predicted one. 

In the last section we interpret the meaning of the approximation used and the discrep- 
ancies that we noticed with the simulations, and point out the direction of future work. 
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2 Closing probabilities 



In the following, we will use the overlap q(C, C) = 1 — d{C, C) to describe the metric of 
configuration space. 

We will follow a stochastic strategy, with the aim to compute the probability distribution 
of the overlap q(t, t') between the configurations C(t) and C(t'), t' > t, on the same trajectory. 
It is convenient to consider only trajectories not yet closed at time t' — 1, so what we actually 
want to compute is a conditional probability, that we denote by the symbol P/v(t, q). 

We define the closing probability ir N (t,t') as the probability to find q(t,t') = 1, which 
means that the two configurations are equal, over a trajectory not yet closed at time t' — 1: 

n N (t, t') = P{q (t, t') = l\ q(h,t 2 ) < 1, h < t 2 < t'} . (3) 

We will have also to deal with the integral closing probability, 7fjv(£), which is the prob- 
ability that a trajectory not yet closed closes at time t regardless of the length of the cycle 
reached, 7f jv(^) = J2i 7Tjv(i — I, t). 

Then it is easy to compute Fjy(t), the probability that a trajectory has not yet closed at 
time t. It satisfies the equation F^{t + 1) = -F/v(t) (1 — Tr^it)), whence, in the limit of large 
systems, introducing a continuous time variable, and using the fact that 7fjv(t) is small even 
at times when nearly every trajectory has closed, one finds the simple expression 

F N (t) = exp (- J* n N (r)dr^j , (4) 

so that the probability to find a trajectory that, after a transient time t, reaches a cycle of 
length /, is given by 

P {T = t, L = 1} = exp ( - f +l n N (r)dT ) n N (t, t + l). (5) 



3 Annealed approximation: the overlap distribution 

The goal to compute the overlap distribution becomes easier if one can reduce the dynamics 
of the overlap to a stochastic process. 

We can treat Kauflman model as a stochastic process by extracting the values of the 
response functions not at the beginning, but every time a variable receives a signal that 
it had never received, and using the values already extracted every time the same signal 
comes back. Note that, if every gene receives signals from the whole system, we have not 
to worry about this last eventuality, which never occurs before the trajectory closes. In this 
case the annealed approximation that we are going to define is exact. This case corresponds 
to the limit of infinite K and is known as the Random Map model, studied ananlitically by 
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Derrida and Flyvbjerg ||. In this work, we generalize their results to finite K values, in the 
framework of the annealed approximation. 

Thus, in order to compute the overlap between two configurations, for every gene we 
distinguish two cases: either the gene received the same signal Si(t) at the two previous time 
steps, and so it will be in the same state, or it didn't: 

q(t + 1, if + 1) = — + i 1 ~ S Si(t)Si(t')) S m(t))f(Si(t'))- (6) 

i i 

Now, the annealed approximation consists in this: we extract at random at every time 
step all the connections in the network, keeping memory only of the minimal information 
about the overlap q(t,t'), so that the probability that the two signals Si(t) and Si(t') are 
the same is q(t,t') K ; if they are not equal, we extract at random the values of the response 
function, and they will be equal with probability 1 — p, where p = 2p(l — p) @]. 

Derrida and Pomeau introduced the annealed approximation in order to obtain the equa- 
tions for the evolution of the average Hamming distance between two randomly chosen con- 
figurations; we will use it in order to treat the overlap as a Markovian stochastic process, 
and so to compute the closing probabiliies. We will argue later from our simulations that 
the overlap distribution obtained through the annealed approximation is very close to the 
quenched one. 

It is not possible, using this kind of approximation, to impose the condiction that the 
trajectory where we measure q(t,t') was not closed at time t' — 1: we can only impose that 
it did not close on a cycle of length I — if — i. As we shall see, this fact is likely to be the 
source of the main discrepancies between the annealed approximation and our simulations. 

If we impose this condiction, we obtain the master equation 



P N (t + l,t + l + l ]qn )=^ n j l-P N (t,t + l;l) ' (7) 

where q m = m/N and 

1 {x) = {l-p)+px K . (8) 

Given an initial distribution Pjv(0, l',q), the process defined by (0), which is ergodic 
because we exclude the overlap q = 1, evolves to a stationary distribution independent 
on both t and I, which appears only in the initial distribution but not in the transition 
probability. 

We look for a stationary distribution in an exponential form, 

P N (q = x) = C N (x)e- Na ^. (9) 
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In the limit N — > oo we can use continuous variables, transform the sum into an integral 
and apply the saddle point method to perform it, obtaining the following selfconsistence 
equations for the unknown functions a(x) and q(x) (the last one is the saddle point of the 
integral for a given value of x): 

a'(q(x)) - i{q{x)) - = (10) 

a(x) = a(q(x)) +x\og (^y) + (1 - x) log ( jz^j ) , (11) 

with the condition q(x) < 1. 

Deriving the first equation and inserting the result in the second one we get a non-local 
transcendent equation for the function q(x): 



if , \\{ x l ~ x \ i ( l( x ) \ i I 7 (<? (l( x ))) \ / 10 x 

We were not able to solve this equation, so we had to solve numerically a discrete version 
of it in order to obtain the exponent of the closing probability, a(l). 

Before to turn to this point, however, there are some important features of the overlap 
distribution that can be easily understood. 

Deriving (|TT|), we can see that the point Q(t) where the overlap distribution is concen- 
trated satisfies the equation 

Q(t+1)= 7 (QW), (13) 

which, at stationarity, becomes 

Q* = 1 (Q*) = l-p + pQ* K . (14) 

This equation was first derived by Derrida and Pomeau as a mean field solution of the 
annealed model [|J. The condition 7'(1) = Kp = 1 separates parameter space in two different 
regions: 

• Kp < 1: frozen phase. The only fixed point of the map (|13D is Q* = 1. This implies, 
as we shall see, that a(l) goes to zero and the typical cycle length increases less than 
exponentially with system size. 

• Kp > 1: chaotic phase. The fixed point Q* — 1 becomes unstable (moreover, it would 
give raise to a negative variance) and the overlap is concentrated around another fixed 
point of dl3|) , Q* < 1. In this case a(l) is finite and the typical cycle length increases 
exponentially with system size. 
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The variance of the distribution can also be easily obtained, in the Gaussian approxi- 
mation, deriving flllD twice and imposing the saddle point condition. One finds that the 
variance multiplied by N, V(t), follows the recursive equation 

V(t + 1) = (V (Q(t)f V(t) + Q(t + 1)(1 - Q(t + 1)), (15) 
whose fixed point is 

v , = Q'(i-Q-) (16) 

It is appearent that the stationary overlap distribution is a non binomial one: the presence 
of correlations between the genes are made evident by the fact that V* is larger than Q*(l — 
Q*). Moreover, it is also easy to see that a'(l) is approximately equal to pK and does not 
diverge, while in the binomial distribution it does. 

On the other hand, when K tends to infinity destroing the correlations between con- 
secutive overlaps, the stationary overlap distribution becomes a binomial one and Q* tends 
exponentially to 1 — p. 

In the frozen phase, when 7'(1) < 1, Q* = 1 is the true solution of equation (|HD, and 
a(l) vanishes or, to say better, is of order 1/N. 

It is important to see how the term of 0(1) vanishes with time: to see that, let us see 
what happens at the critical point p = 1/K. 

The evolution of the most probable overlap at the critical point, for K > 1, is such that, 
introducing a continuous time variable, 

Q'(t) = hQ-l/j:(K-l-z)Q\ (17) 
^ i-i 

and the asymptotic behaviour in t is 

i-Qocf^r. (is, 



Let's now turn to the numerical computation of the exponent of the closing probability, 
«(!)• 

We have to use the discrete variables = i/N; moreover, in place of the stationary state 
equation, we consider the evolution of the exponent of the overlap distribution, a t ({qi}) (we 
omit to specify the dependence on I, which appears only in the initial condition ao({qi}))- 

Neglecting the normalization factor, in the master equation, we get the following map 
for the evolution of a t (qi): 

ait+i (qi) = qi log (ft) + (1 - ft) log (1 - ft) + 

min j<N (a t (qj) - (ft) log (7(ft)) - (1 - ft) log (1 - 7(ft))) , (19) 



7 



This map does not have a fixed point {a(qi)}, as it can be seen numerically: for istance, 
cxt(l) reaches a minimum at a certain time and then increases uniformly. 

This increase is not in contradiction with the existence of a stationary distribution; in fact 
equation flTUP misses the normalization factor (1 — P(l))~ , and so the probability computed 
by means of it has to decrease when the exact probability has reached the stationary value. 

The increment is of order O (^J, because, after Qsit) reaches its stationary value Q*, 

a(Q* N ) grows of a term of that order due to the difference, O (^j , between Q* N and 7(Qjv), so 

it is not important because we also miss terms of order in the saddle point approximation. 
In our algorithm, we decided to get the function a t (qi) at the time when a t (l) reaches its 
minimum value. 

There is another finite N effect in our numerical calculation: because of the condition 
j < N (this means that the overlap Xj = 1 cannot be taken as a starting point), the 
stationary value of a^ixi) decreases with N. In fact, the larger is N, the lower can be the 
minimum in (|T9|), and the effect, of order 1/N, is important especially for i = N, because in 
this case in the stationary state the minimum in ( |I9| ) would be attained for j = N. 

The finite size effects become more relevant at the critical line p — 1/K. As a discretiza- 
tion effect, the most likely value of the annealed distribution is not Q* = 1, because the 



corresponding value j = N in the equation (|T9|) is not allowed, and at stationarity a{xi) has 



a minimum for Q* N = 1 — e N , where is O (ttW)j so that i(Q*n) ~~ Q*n * s °f or der O (^)- 
Actually, in the frozen phase the equation ( |19|) is not more a good approximation, because 
then we are not allowed to neglect -P(l) in the denominator of equation (0). We keep it in 
the complete form and obtain the equation for the first moments of the distribution: 

(I - p) + p(q K ) = (q) (I - P(l)) + (P(l)) , (20) 
whence, at the critical point, neglecting terms of order ((q — (q)) 3 ) and (1 — (q)) 3 , one gets 

(21) 

This suggests that both -P(l) and 1 — (q) are of order 1/y/N. Combining this fact with 
the previous results, one finds in the critical line an asymptotic closing probability of the 
form 

/ i\ 1 / -4iV \ . . 

n N (t, t') cc -= exp — — -— , (22) 




while in the chaotic phase we find an exponential behaviour: 

7t N {t,t') <xe a ^ N , (23) 

provided that the smallest time t is large enough for the overlap to reach the stationary 
distribution. 
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4 Period distribution within the Annealed Approxima- 
tion 

We now use the closing probabilities obtained in the previous section to compute the asymp- 
totic distribution of transient times T and cycle lengths L. 

In the chaotic phase, inserting the expression for the asymptotic closing probability into 
the equation (|[), we obtain 

Prob {T = t, L = 1} ~ exp + l) 2 e~ aN ^j e~ aN , (24) 

(from now on, we will shortcut a(l) with a). 

So, for finite K, we obtain the same form of the distribution as in the Random Map 
Model, differing only in the value of a. 

Introducing the continuous variables x = t/r and y = l/r, with 

r = e^ aN , (25) 
the period and transient distribution becomes a continuous distribution with density 

f(x,y) = e-^ 2 . (26) 

While the form of the distribution predicted by the annealed approximation is not con- 
firmed by our simulations, the N behaviour of the time scale r is in good agreement with 



the results of the simulations, as it is shown by figure [12] where we compare the value of 
a obtained by the iteration of equation (|19|) with twice the exponent of the average cycle 
length and transient time or with the exponent of the variance of the period obtained by a 
fit of experimental data, for some K value. Other numerical results will be presented later. 



In the frozen phase a vanishes, and the average cycle length grows as a power of N. The 
closing probability, for K ^ 1, is given by 

ir N (t, t + I) oc ^. ]_ exp (- ) ^, (27) 

where the time-scale at the critical point is given by 

r = JN/(K-1). (28) 

This suggests a kind of universality of the critical line in the Kauffman model: the time 
scale of the dynamics diverges always as the square root of the number of the elements 
(except for the case K = 1, where the variance in the Gaussian approximation diverges 
linearly with t, thus destroying the validity of the approximation, and our method fails). 
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It is easy to see that the exponent a of the closing probability vanishes quadratically when 
p approaches to the critical point (except for the case K — 1). In fact, putting p = p c + e, 
one finds 

IK 

3* = 1-^33-6 + o(e), (29) 



1 (i - Q*f my 2 

a ~ ~ = — — -e 



(e 2 ) (30) 
e of infinite K, where the o 

distribution is a binomial one with mean value 1 — p and a is equal to — log(l — p)). 



2 V K-l 

(of course these expressions are not valid in the limit case of infinite K, where the overlap 



5 Distribution of weights 

In Kauffman networks configuration space breaks into a number of limit cycles with their 
own attraction basins. The statistical distribution of the weights of such basins, W a (the 
rate of the number of configurations in the basin to the total number), has been computed 
in the limit cases K — 1 and K = 00 by Derrida and Flyvbjerg [[J and Flyvbjerg and Kjaer 
fTTH respectively 

In the framework of the annealed approximation it is possible to obtain this distribution 
for every value of the parameters, at least in the chaotic phase: in this case, it turns out to 
be the same distribution computed in || for the Random Map. 

Following 0, it is convenient to start computing the quantities 

% = YMi (3i) 

a 

that are the moments (of order p — 1) of the probability to find a basin of weight W. They 
can be computed by noting that they represent the probability that p configurations chosen 
at random ultimately meet. 

In the framework of the annealed approximation, one finds that these quantities do not 
depend, in the chaotic phase and in the limit iV — > 00, on the parameter r which gives the 
cycle length time scale, and have the same values computed by Derrida and Flyvbjerg for 
the case of the Random Map. 

To be more concrete, let's calculate Y2. We must generate randomly two trajectories 
on the same network. We extract the first one and follow it till the time Ti when it closes 
on itself, as we have done in the previous section; then we generate another "replic" of the 
system and we study the overlap ^12(^1,^2) between the configuration at time t\ on the first 
trajectory and the one at time £2 on the second one. 

We define the closing probability of the second trajectory on the first one as the proba- 
bility that ^12(^1,^2) is equal to 1, with the condition that the first trajectory is composed of 
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Ti different configurations, the second of t 2 — 1 different ones and none of them has met the 
first trajectory before time £ 2 (we express this condition by the symbol 0(T\, t 2 )): 

■KfiTiMM) = P{qi2{txM) = l\0{T x ,t 2 )}. (32) 
The probability that the second trajectory meets the first one after T 2 time steps is then 

4 2) (T l5 T 2 )4 2) (T i; T 2 ), (33) 

where 7fjy (Ti;i 2 ) is the integral closing probability, say the sum over t\ of 7% (Ti;ti,t 2 ), 
and Fjy (Ti,T 2 ) is the probability that the second trajectory is composed of T 2 different 
configurations and they don't touch the first trajectory before time T 2 . One has, in the limit 
of continuous time variables, 

F§\T U T 2 ) = exp (- J^' 1 (^ 2) (T i; t) + 4 2) (t)) dt^j , (34) 

and 

Y= F N (T 1 )F$\T 2 ,T 1 )n N (T 1 )n% ) (T 2 -,T 1 ), (35) 

Ti,T 2 

where 7fjv(i) an d T\r(t) are respectively the integral closing probability and the opening 
probability for only one trajectory, as defined in the second section. 

In the framework of the annealed approximation, we can imagine that the overlap qi 2 (t, 0) 
between the initial configuration of the second trajectory and a given configuration of the 
first one is a binomial variable, and that the distribution of q\ 2 {t + 1, 1) is obtained from it 
through the annealed evolution equation. 

In this hypothesis, asymptotically in t\ and t 2 , the closing probability in the chaotic 
phase is given by 7rjy (Ti; ti, t 2 ) = 1/t 2 , where r = exp(aiV/2) is the cycle length time scale. 

Operating the substitutions x = Ti/r.y = T 2 /t and transforming to continuous variables, 
one then finds, for large values of r, 

Y 2 ~ / x 2 e —dxdy = 2/3, (36) 

Jo Jo 

which is independent on r and is the same result of JJ. 

This computation generalizes easily to an arbitrary number of trajectories, and in this 
way one finds for everyone of the Y p the same value obtained by Derrida and Flyvbjerg for 
the Random Map; so one can argue that f(W), which has the meaning that f(W)dW is 
the average number of cycles with weight between W and W + dW, is given, in the chaotic 
phase and in the framework of the annealed approximation, by the same expression valid in 
the Random Map model: 
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f(W) = -W- 1 (l-W)- 1/2 . (37) 

To compare this prediction with numerical data, we have pursued the simulations per- 
formed by Derrida and Flyvbjerg || relatively to Y 2 for K = 3, doubling the size of the 
maximum system simulated. For each system size we generated 10000 sample networks and 
two configurations at random on each of them and we measured Y 2 as the probability that 
they end up on the same cycle, as it was done in ||. 

In this way we obtained a curve of Y versus N which reaches a minimum value of 
Y = 0.590 at N = 40 and then increases with N. The largest system that we could simulate 
(N = 100) looks to be still far from having attained the stationary value, so we are unable 
to say if it agrees with the annealed prediction Y = 2/3, but it is very close to that value 
(figured). 

We measured also the histogram of the probability to find a cycle of weigth W. Figure 
|2| compares the histogram obtained generating at random 1000 networks with K = 3 and 
N = 50 and simulating 250 trajectories on each of them with the one computed using (|37|), 
which is valid in the infinite N limit. The agreement, as it is possible to see, is not so bad. 



6 Average cycle number 

In this section we aim to compute the average number of cycles in a network, a quantity 
introduced by Kauffman to represent, in his biological metaphor, the typical number of 
possible cell types for a genome of given length. This quantity is difficult to observe in 
the simulations, because one must consider also cycles of very small weight, and has been 
computed only for the Random Map (infinite K, p — 1/2; ||). 

For this purpose, we have to compute the distribution of the overlap between a configu- 
ration at time t and the initial configuration, P/v(0,t; q), which plays the role of the initial 
distribution for the stochastic process defined by the master equation (0). 

It is easy to generalize the annealed approximation in order to get the transition prob- 
ability of the joint distribution of the overlaps q{0,t) and q(t — l,i): every gene remains at 
time t + 1 in the same state it was at the previous step with probability 7 (q(t — 1, t)), where 
7(x) is defined in the previous section (||). 

If the genes unchanged are Nz of the genes whose state was the same in the initial 
configuration and N(q(t,t + 1) — z) of the other genes, the overlap q(0,t + 1) is given by 
1 — g(0, t) — q(t, t + 1) + 2z. One then obtains the transition probability 

P {q(0, t+l)= Xq, q(t, t+l)=Xi\ q(0, t) = q , q(t - 1, t) = qx} = 

(IS) i^)f X1 (1 - ^)) N(1 ~ X1) , (38) 
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Figure 1: Y 2 versus N in networks with K = 3. The annealed prediction is Y 2 = 0.6. 
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Figure 2: Histogram of the probability to find a cycle of weigth W in networks with K = 3 
and N = 50, obtained simulating 1000 networks and 250 trajectories on each of them. The 
dotted line is the same thing, computed in the framework of the annealed approximation. 
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with 

z = - (xi + x + g - 1) , (39) 

and with the bounds 

|1 - (x + xi)\ < q < 1 - \x -xi\ . (40) 

At time t — 1, the variables q(0, t) and q(t — l,t) coincide and their distribution is a 
binomial one, but the correlations disappear after 2 time steps, at least in the Gaussian 
approximation, and the stationary joint distribution factorizes into the stationary distribu- 
tion of q(t — l,t), that we met in the previous section, and a binomial distribution of mean 
value 1/2. The most likely value of q(0,t) is always 1/2, even before that the stationary 
distribution has been reached, and the exponential part of the closing probability 7rjv(0,t) 
(the probability to have q(0,t) = 1) is always 1/2 N , but the prefactor is different from 1, in 
the initial steps. As usual, this can be seen putting the distribution in the exponential form 
and performing the integral with the saddle point method. 

The average number of cycles of length / is related to the probability P(l) to find a 
configuration belonging to a cycle of length / by the obvious formula 

fn = ^P{L = l,T = 0}. (41) 

In the framework of the annealed approximation, using the results of section 4 about 
period and transient distribution and the fact that 7^(0,/) = 1/2 n q one finds, in the 
chaotic phase, 

Q / U 2 ' 



n, = -iexp(--— . (42) 



Summing over / one gets the average number of cycles in a network. Apart from systems 
at the critical point, c\ is a quantity of order 1, and tends to 1 when / grows; so, substituting q 
by 1 and transforming the sum into an integral, one finds, asymptotically in r, the following 
behaviour: 

2^ 1 

Y,n t w logT = -aN, (43) 
i=i 1 

where a is the exponent of the stationary closing probability introduced in the previous 
section. 

This formula is valid only in the chaotic phase: in the frozen phase, where a cancels, we 
are not allowed to put q to 1. 

For instance, for I — 2 one finds 

MO, 2) = ^f e-^^l-i))" 1 ^, (44) 
2 J J2n/N 
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where 



f(x) = xlog(x) + (1 — x) log(l — x) — x\og{ n f{x)) — (1 — x) log(l — 7(2;)) ■ (45) 

(7(2;) is defined in (|8|)). The function /(x) is concentrated around the same point Q* around 
which it is concentrated the stationary distribution of the overlap, i.e. the solution of the 
equation Q* = j(Q*). The saddle point method with Gaussian corrections yields the result 

M0,2) = ^(l-7'(Q*)r\ (46) 

but at the critical point j'(Q*) is equal to 1, so one has to go beyond Gaussian corrections 
and finds 

MO^oc^A^. (47) 

7 Numerical tests of the Annealed Approximation 

All simulations reported in this section are referred to Kauffman Networks with unbiased 
response functions, i.e. the parameter p is fixed at the value 0.5, and we will no longer 
mention it. 

Our numerical analysis of the overlap distribution confirms the validity of the annealed 
approximation in some range of the time variables and reveals deviations from it in some 
other range, so that the period distribution turns out to be different from the one predicted 
by the annealed approximation, while the time-scale of the period and the first moment of 
the distribution of the weights are in good agreement with our computation. 

7.1 Overlap distribution 

We measured the overlap q(t,t + /) on the subset of the trajectories not yet closed at time 
t + 1 — 1, in order to sample with the conditional probability defined in section 2. In this 
way we also eliminate the "deterministic noise" due to the periodicity of the orbits. 

For each value of the parameters, we generated 20 '000 sample networks and on each of 
them we simulated a trajectory measuring the overlap q (C(t), C(t + I)) till the time when 
the trajectory closed. So we obtained the histogram of the overlap q(t,t+l) for chosen values 
of the temporal distance / and t large enough to suppose that the distribution is stationary. 

From the annealed approximation we expect the stationary distribution not to depend 
on I, but in fact we could observe different behaviours for different I values. So the data 
analysis is rather complex since for every point (K, p) in parameter space three varables, N, 
t and /, are involved and different limits, e.g. N large at fixed t and / or t and I large at 
fixed N, may give different behaviours. 
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We are interested to the exponential part of the stationary distribution, and we estimate 
it through the formula 



oci{q) « ^ {--logN - \ogP m {q)^ , (48) 

in which the factor 1 / coming from the Stirling expansion of the binomial coefficient has 
been taken into account (actually, we use this formula only for 0.1 < q < 0.9, while for q = 
and 1 we don't subtract the term with logiV and we interpolate linearly for intermediate 
values of q). 

The function of q so obtained is compared with the function a(q) numerically computed 
in the annealed approximation (equation [19]). 

Most of our simulations have been performed on systems with connectivity K = 3. 

The agreement is very good when I, the temporal distance between the configurations, 
is large (in fig. [3] (b) data are shown for I = 62, in networks with K = 3 and iV = 50), while 
the shape of the distribution is different, especially close to q — 1, for small values of I (in 
fig. H] (a) we show data for / = 2, in networks with K = 3 and N = 50). 

In the frozen phase the agreement between the annealed approximation and our data 
worsen. We simulated systems with K = 2 and different values of N. We show in figure 
01 the function a(x) computed numerically within the annealed approximation and the one 



measured in the simulation with N = 90 using formula (48) (in this case, we subtracted the 
log term also for q = 1, because we expect P(l) to be proportional to 1/yN). 

The agreement between data and computation is again better when the temporal dis- 
tance between the configurations compared is large, but the real distribution is always much 
broader than the annealed one and the small overlap values are much more likely in the 
simulation than in the annealed approximation. 



Then we tried to see if the overlap between configurations on the limit cycles is sta- 
tistically different from the overlap between transient configurations. For such purpose we 
computed the function ai(q) imposing two kind of conditions: on one hand, that both con- 
figurations be on the cycle, on the other that both configurations be transient. 

In the chaotic phase (we report data with K = 3 and N = 50) the functions so obtained fit 
within one or two standard deviations (which are large, because the selection of trajectories 
is very severe) to the annealed prediction, for the same values of parameters for which the 
overall distribution gives a good function a(q), and so we don't see differences in the overlap 
distribution between pairs of transient configurations and pairs of configurations in the limit 
cycles. 

The situation is different in the frozen phase where one can see that the most probable 
value of the overlap between transient configurations is quite lower than the most probable 
overlap between configurations of a cycle. 
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Figure 3: The function ai(q), which is the exponential part of the distribution of the overlap 
q(t,t + I). The solid line is the annealed prediction, the points plotted come from the simu- 
lation o/20000 networks with K = 3 and N = 50. The distribution is obtained considering 
different values oft between 50 and 80. (a): I = 2; (b): I = 62 (in this case the agreement 
for q close to 1 is much better). The condiction that the trajectory must not be closed at 
time t + / — 1 selects respectively 79 an^ 63 per cent of trajectories in the sample. 




Figure 4: The function cti(q), which is the exponential part of the distribution of the overlap 
q(t, t + l). The solid line is the annealed prediction for K = 2 and N = 90, the points plotted 
come from the simulation of 40000 networks with K = 2, and N respectively equal to 60 
(diamonds) and 90 (asterisks). The distribution is obtained considering different values oft 
between 12 and 20. The vertical scale is logarithmic, (a): I = 2; (b): I = 26. The condiction 
that the trajectory must not be closed^ time t + 1 — 1 selects respectively 59 and 12.5 per 
cent of trajectories in the sample for N = 60 and 76 and 19 per cent for N = 90. 




Figure 5: The function cn(q) again. Here K = 2 and N = 90. The solid line is the annealed 
prediction, the points plotted come from the simulation of 40000 networks. The distribution 
of the overlap q(t,t + l) is obtained considering different values oft between 12 and 20. The 
vertical scale is logarithmic. The two graphs concern 1 = 2 (a) and I = 14 (b) respectively. 
The asterisks are refered to the overlap of two configurations in a limit cycle, the diamonds 
to the overlap between transient configurations. 




Figure 6: The function cn(q) again. Here K = 3 and N = 50. The solid line is the annealed 
prediction, the points plotted come from the simulation of 20000 networks. The distribution 
of the overlap q(t,t + l) is obtained considering different values oft between 50 and 80. The 
two figures are concerned with 1 = 2 (a) and I = 62 (b) respectively. The asterisks are refered 
to the overlap of two configurations in a limit cycle, the diamonds to the overlap between 
transient configurations. 91 



7.2 Measure of the closing probabilities 



In section 2 we denned the closing probability 7rjy(t, t + I) as the probability that the con- 
figurations at times t and t + I be equal, with the condition that the trajectory is not yet 
closed at the previous time steps. 

We studied the closing probability for I fixed as a function of t, and we saw that in 
the quenched model the overlap distribution doesn't go to a stationary distribution as it is 
predicted by the annealed approximation. 

For a fixed system size, the closing probability grows to a maximum value, usually larger 
than the annealed stationary value (but this depends on /!) and after it is always decreasing. 

The rate of decrease appears to be independent on I, while the maximum value attained 
is not: it is larger if I is even and decreases with / for a given parity. 

In figure [7| we report data for K = 2, which, in our case p = 0.5, is the starting point of 
the frozen phase, and K = 3 in the chaotic phase. 

The simulation results are compared to the annealed closing probability exactly computed 
iterating equation (|^). The time behaviour is at first in good agreement (especially for I 
large), over a fairly large range: in the networks with K = 2 and iV = 120 the annealed and 
the quenched grow together from 10 -36 to 10~ 2 . After, the annealed remains constant while 
the quenched decreases. 

On the other hand, when we compared to the annealed approximation the closing prob- 
ability measured at fixed times varying system size, we found out that the agreement, at 
the beginning not very good (many standard deviations), becomes very satisfactory taking 
larger systems. 

Put into a graphic, the N behaviour of the closing probability appears to be an expo- 
nential decay with small corrections. 

We now study the integral closing probability, Tf N (t), which is the sum over t' of 7rjv(i', t). 
As usual, we considered two values of K, 2 and 3, in and near to the frozen phase, and iV 
ranging from 60 to 180 for K = 2 and from 20 to 50 for K = 3. 

Also vfAr(t) has a non monotone behaviour for t not too large: it increases to a maximum 
value and then start decreasing. Of course, it has to increase at large times, because it has 
to fulfil the normalization condition tt (% n ^ = 1, but in our simulations it keeps decreasing 
at all times accessible to measurement. In the framework of the annealed approximation, we 
expect that it grows linearly with t. 

There is nothing new in this behaviour, which is just a consequence of the decrease in time 
of the closing probability, but it is the origin of other interesting features of the quenched 
system. In a next subsection we shall try to give an interpretation to this funny result, which 
is in contrast with the annealed approximation. 

We show the temporal behaviour of 7Tat(£) for different system sizes in the chaotic phase 
(fig. ^) and in the frozen one (fig. |8|). The two plots have qualitatively the same shape, but 
in the frozen phase the integral closing probability decreases very slowly with system size 
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Figure 7: Closing probability 7T^(t,t + I). In solid it is reported the annealed value with 
1 = 1. The ordinate scale is logarithmic, (a): The system is in the chaotic phase, K = 3 and 
N = 36. Data correspond to I = 1 (middle), I = 2 (upper line: the closing probability is sys- 
tematically higher) and to I = 15 (diamonds; the agreement with the annealed approximation 
is better), (b): The system is at the critical point, K = 2 and N = 120. Data correspond 
I = 2 (upper line) and to I = 15 (diamonds). The first point of the annealed computation, 
P = 7.52 • 10~ 37 ; is not shown. 




Figure 8: Integral closing probability, vrjv(t) in the frozen phase (K = 2) for different system 
sizes: N = 60, 90, 120, 150, 180. n N reduces when N grows. 



24 




Figure 9: Integral closing probability, 7Tjv(0 in the chaotic phase (K — 3) for different system 
sizes: N = 20,30,40,50,60. tt n reduces when N grows. 
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N, while in the chaotic phase it decays exponentially. 

We also measured the mean overlap between successive configurations, q(t, t+1), without 
limiting the sample to the only trajectories not yet closed. We considered three different 
connctivities in the chaotic phase (K=3,4 and 5) and varied N. 

As a function of t, this quantity increases to a stationary value which tends, when N 
grows, to the same value Q* computed within the annealed approximation, but the (positive) 
corrections are large even for system sizes of the order of hundreds; on the other hand, when 
the temporal distance between the configurations compared is large, the corrections are much 
smaller. 

This measurement was done without imposing the opening condition. When we mea- 
sured the mean overlap between successive configurations, q(t,t + 1), on the subset of the 
trajectories not yet closed at time t + 1, we noticed that this quantity has a non monotone 
behaviour. It starts at t = with the value 0.5, then increases, reaches a maximum value 
at about the same time at which the integral closing probability has its maximum value, 
and decreases again. Thus, the behaviour of the mean overlap parallels that of the closing 
probability. 

This behaviour indicates that the condition that the trajectory is not yet closed at a 
certain time selects, from a point on, smaller and smaller overlaps. Fig. |ll| shows the mean 
and the variance of the distance d = 1 — q of successors on trajectories not yet closed and 
the integral closing probability for K = 2 and N = 90. The effect is more evident in the 
frozen phase, but it is still present in the chaotic phase for low values of K. 

7.3 Properties of the cycles 

Our simulations confirm well the predictions of the annealed approximation concerning the 
time scale of the period and of the transient distribution. 

In the chaotic phase these quantities have an exponential behaviour, exp(aiV/2). We 
simulated systems with K ranging from 3 to 6 and several values of N (from 10 to 100 for 
K = 3, from 16 to 40 for K = 4, from 16 to 31 for K = 5 and from 6 to 28 for K = 6), 
generating at least 10000 samples of each network and few configurations on each network 
(except for the case K = 3, where we generated hundreds of initial configurations in order to 
compare the variance of cycle length on a given network with the variance between different 
networks) . 

We fixed a maximum time for the simulation, exponentially increasing with N; in such a 
way, we underestimated the average cycle length because some trajectories where left before 
they completed their limit cycle; they were about one percent in the worst cases, but we 
think that this does not affect very much the determination of the exponents of the average 
period. 

In fig.|12| we plot the exponents a of the time scale of the system, defined by r = 
exp (aN/2), and compare it with the annealed prediction, numerically computed. 
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Figure 10: (a): Average overlap between successive configurations The solid line shows the 
annealed equation, Q(t + 1) = ^(l + Q(t) K ^j. (b) Variance of the overlap between suc- 
cessive configurations. The solid line shows the annealed equation, V{t + 1) = Q(t + 
1) (1 — Q{t + 1)) + {^Q(t) K ~ l ^j V(t). The system is in the chaotic phase, K = 3; data 
are shown for N = 30, 80, 200 and 400; the larger is N, the lower is the curve. 
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Figure 11: Average ((}) and variance (+) of the distance d = 1 — q between successive 
configurations on trajectories not yet closed at time t — 1, as a function oft. It is also shown 
the integral closing probability, 7Tat(0 (•)■ K = 2, N = 90, 10000 sample networks. 
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It is possible to measure in the simulations four different exponents: ai (twice the mean 
period exponent), ay (the exponent of the variance of the period), olt (twice the mean tran- 
sient exponent), and ap, the exponent of the time-scale of the period distribution (defined 
by fitting the probability to find a period larger than t with the function exp ( — (t/r) 13 )). 



These exponents were obtained through a fit of data with different iV values. Their values 
lie very close whithin a few percent. In the figure there are not error bars, because the fits 
of our data are affected by finite size effects whose importance is difficult to estimate, but 
the agreement is quite good, even with not very large systems. 

The largest deviation was found for the case K = 3, where the annealed approximation 
gives a = 0.182 while our simulations give ol = 0.210, ax = 0.204, ap = 0.21 and ay = 
0.206 for the average variance in a given network and 0.205 for the variance between different 
networks. 

In this case, we performed the annealed computation with the map ( |I9|) and the same 
values of N used in the simulations, in order to obtain the corrections due to the discretiza- 
tion, whose effect is to increase a. In figure |H| we compare twice the logarithm of the average 
period (diamonds) with Na(N), where a(N) is computed with the annealed approximation. 
The discrepancies now are reduced respect to the asymptotic a value, but the simulation 
data still seem to be systematically larger and the angular coefficients are different. 

In the simulations with K = 3 we measured also F^(t), defined in the second section 
as the probability that a random chosen trajectory is not yet closed at time t, and related 
to the closing probability by formula (|j). According to the annealed approximation, one 
should find log(F/v(t)) oc t 2 , because the closing probability reaches a constant value. In 
fact one finds that this quantity is linear or less than linear in t, so that Fjsr(t) has a stretched 
exponential behaviour. 

This feature is not surprising because we have already seen that the integral closing 
probabiity is a decreasing function of t from a point on, and log (F/v(t)) is just its integral, 
but it is striking that the exponent of t is a slightly decreasing function of N. 

This last observation was made looking at the period distribution which is a quantity 
much easier to compute: to measure F^lt) one has to keep in memory the whole trajectory, 
while to measure the period one needs only a reference configuration. 

The distribution of cycle length L and transient time T are simply related to -Fjv(t): 



For K = 3, the probability to find a cycle of period greater than t decays as a stretched 
exponential, exp (— (t/r)^j; the time scale r grows exponentially with N, with the same 
exponent a^ of the average cycle length; the decay exponent (3 goes down with N, even if 
very slowly: we found an exponent of 0.91 for N = 10 and 0.39 for iV = 100. This fact is in 
contradiction with the results of the annealed approximation, which predicts (3 = 2, but it 
is a direct consequence of the decrasing with time of the closing probability. 



Prob {T = t,L 



1} = F N (t + l)% N (t,t + l). 



(49) 
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Figure 12: The function ot{K), where a is twice the exponent of the characteristic time-scale 
of period and transient distribution, r = exp (aN/2) , for p = 0.5, as obtained from the 
annealed approximation. The asterisks are obtained from a fit of simulational data; the last 
one is the theoretical result for the Random Map. 
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Figure 13: ((}) : twice the logarithm of the average period in networks with K = 3 and N 
nodes (simulation); (+) : the annealed coefficient a computed in the saddle point approxi- 
mation for a system of N nodes, times N . 
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From our simulations on systems with K = 2 and iV ranging from 30 to 90 and K = 3 
and N ranging from 10 to 100 it appears that the histogram of transient time has a peak 
corresponding to the time t where the closing probability reaches its maximum value, and 
then decreases as a stretched exponential; and the cycle length histogram decreases with the 
same law. 

In addition, very short periods are much more common respect to the others than one 
would expect on the basis of the annealed approximation (see figure 14). 

Another unexpected feature of the period distribution is the fact that cycles of even 
lengths are more likely than odd ones, as one can already see by the fact that the closing 
probability is higher for even periods. 

We have seen in this fact a clue of a non trivial distribution of local periodes, say the 
period of a single element in a network. We have obtained in our simulations an indirect 
confirmation of this guess. We will present these data in a forecoming paper. 



8 Discussion 

In this work we used a stochastic scheme, generalizing the analysis by Derrida and Flyvbjerg 
of the Random Map Model to finite values of K, in order to compute the distributions of 
cycle lengths, transient times and attraction basins weigths in Kauffman Networks. The 
definition of the closing probabilities is crucial in this scheme. 

The annealed approximation, introduced by Derrida and Pomeau, has revealed to be a 
good tool for the computation of the closing probabilities and the approximate solution of 
the model. 

There are deviations from the approximation, which are larger close to the frozen phase, 
but the mean quantities that we computed are in good agreement with our numerical results. 
The annealed approximation works better in the chaotic phase, but it allows also to predict 
an universal behaviour of cycle lengths along the critical line. We have in plan to check this 
point. 

We haven't, till now, tried to justify the approximation used. The annealed approxi- 
mation allows to compute the overlap distribution using only a minimal information about 
the past evolution of the system and neglecting everything else, so we think that a loss of 
memory of the history of the system is needed for its validity. This fact is consistent with our 
numerical results, which show that the overlap of two configurations very far apart along the 
trajectory is in some sense "more annealed" than the overlap of two configurations tempo- 
rally closer. How this loss of memory comes about and what are its limitations would be an 
interesting question to address to understand in a more analytical way dynamical systems 
with quenched disorder. 

As we saw in the previous section, the main discrepancies between our simulations and 
the annealed approximation arise from the fact that the closing probabilities ultimately 
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Figure 14: Histogram of cycle (solid line) and transient (dots) length for a system in the 
frozen phase, K = 2 and N = 120. 10000 sample networks were generated. 
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Figure 15: Histogram of cycle (solid line) and transient (dots) length for a system in the 
chaotic phase, K = 3 and N = 50. 10000 sample networks were generated. 
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decrease in time. 

This fact can be interpreted as a consequence of the definition of the closing probabilities 
as conditional probabilities: the requirement to consider only trajectories not yet closed 
selects, from a time on, trajectories where the typical overlaps are smaller and smaller. 

We did not impose the opening condition in our computation, so we think that the best 
way to improve the annealed approximation would be to take into account the opening of 
the trajectory in the computation. This purpose requires to keep memory of much more 
information about the past history of the trajectory than the simple one that we used in 
doing the annealed computation. We are planning to try to do this in a next work. 
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